Основные задачи, решаемые в страховании
Риски, оцениваемые в автостраховании
Основная идея метода бустинга
Особенности бустинга "деревянных" моделей
Схема решающего дерева для задачи классификации
Основная мера неоднородности при построении деревьев - энтропия Шеннона
,
где p - вероятность состояния, i - состояние системы. Чем выше энтропия, тем больше хаоса.
Прирост информации (изменение энтропии)
Алгоритм построения решающих деревьев
Для задачи регрессии
Градиентный спуск
Что работает лучше, деревья решений или линейные модели?
# Подключение к Google drive
#from google.colab import drive
#drive.mount('/content/drive')
import numpy as np
import pandas as pd
# Загрузим набор данных
#df = pd.read_csv('/content/drive/My Drive/Colab Notebooks/freMPL-R.csv', low_memory=False)
df = pd.read_csv('D:\Cloud\Git\geekbrains-ml-business\\freMPL-R.csv', low_memory=False)
df = df.loc[df.Dataset.isin([5, 6, 7, 8, 9])]
df.drop('Dataset', axis=1, inplace=True)
df.dropna(axis=1, how='all', inplace=True)
df.drop_duplicates(inplace=True)
df.reset_index(drop=True, inplace=True)
В предыдущем уроке мы заметили отрицательную величину убытка для некоторых наблюдений. Заметим, что для всех таких полисов переменная "ClaimInd" принимает только значение 0. Поэтому заменим все соответствующие значения "ClaimAmount" нулями.
NegClaimAmount = df.loc[df.ClaimAmount < 0, ['ClaimAmount','ClaimInd']]
print('Unique values of ClaimInd:', NegClaimAmount.ClaimInd.unique())
NegClaimAmount.head()
df.loc[df.ClaimAmount < 0, 'ClaimAmount'] = 0
Для моделирования частоты убытков сгенерируем показатель как сумму индикатора того, что убыток произошел ("ClaimInd") и количества заявленных убытков по различным видам ущерба за 4 предшествующих года ("ClaimNbResp", "ClaimNbNonResp", "ClaimNbParking", "ClaimNbFireTheft", "ClaimNbWindscreen").
В случаях, если соответствующая величина убытка равняется нулю, сгенерированную частоту также обнулим.
df['ClaimsCount'] = df.ClaimInd + df.ClaimNbResp + df.ClaimNbNonResp + df.ClaimNbParking + df.ClaimNbFireTheft + df.ClaimNbWindscreen
df.loc[df.ClaimAmount == 0, 'ClaimsCount'] = 0
df.drop(["ClaimNbResp", "ClaimNbNonResp", "ClaimNbParking", "ClaimNbFireTheft", "ClaimNbWindscreen"], axis=1, inplace=True)
pd.DataFrame(df.ClaimsCount.value_counts()).rename({'ClaimsCount': 'Policies'}, axis=1)
#conda install plotly
import plotly.express as px
fig = px.scatter(df, x='ClaimsCount', y='ClaimAmount', title='Зависимость между частотой и величиной убытков')
fig.show()
Для моделирования среднего убытка можем рассчитать его как отношение величины убытков к их частоте.
df.loc[df.ClaimsCount > 0, 'AvgClaim'] = df['ClaimAmount']/df['ClaimsCount']
Класс для общих случаев
class InsDataFrame:
''' Load data method '''
def load_pd(self, pd_dataframe):
self._df = pd_dataframe
''' Columns match method '''
def columns_match(self, match_from_to):
self._df.rename(columns=match_from_to, inplace=True)
''' Person data methods '''
# Gender
_gender_dict = {'Male':0, 'Female':1}
def transform_gender(self):
self._df['Gender'] = self._df['Gender'].map(self._gender_dict)
# Age
@staticmethod
def _age(age, age_max):
if pd.isnull(age):
age = None
elif age < 18:
age = None
elif age > age_max:
age = age_max
return age
def transform_age(self, age_max=70):
self._df['driver_minage'] = self._df['driver_minage'].apply(self._age, args=(age_max,))
# Age M/F
@staticmethod
def _age_gender(age_gender):
_age = age_gender[0]
_gender = age_gender[1]
if _gender == 0: #Male
_driver_minage_m = _age
_driver_minage_f = 18
elif _gender == 1: #Female
_driver_minage_m = 18
_driver_minage_f = _age
else:
_driver_minage_m = 18
_driver_minage_f = 18
return [_driver_minage_m, _driver_minage_f]
def transform_age_gender(self):
self._df['driver_minage_m'],self._df['driver_minage_f'] = zip(*self._df[['driver_minage','Gender']].apply(self._age_gender, axis=1).to_frame()[0])
# Experience
@staticmethod
def _exp(exp, exp_max):
if pd.isnull(exp):
exp = None
elif exp < 0:
exp = None
elif exp > exp_max:
exp = exp_max
return exp
def transform_exp(self, exp_max=52):
self._df['driver_minexp'] = self._df['driver_minexp'].apply(self._exp, args=(exp_max,))
''' Other data methods '''
def polynomizer(self, column, n=2):
if column in list(self._df.columns):
for i in range(2,n+1):
self._df[column+'_'+str(i)] = self._df[column]**i
def get_dummies(self, columns):
self._df = pd.get_dummies(self._df, columns=columns)
''' General methods '''
def info(self):
return self._df.info()
def head(self, columns, n=5):
return self._df.head(n)
def len(self):
return len(self._df)
def get_pd(self, columns):
return self._df[columns]
Создаем класс-наследник, в котором переопределяем некоторые методы, специфические для конкретной ситуации, и создаем новые
class InsDataFrame_Fr(InsDataFrame):
# Experience (weeks to years)
@staticmethod
def _exp(exp, exp_max):
if pd.isnull(exp):
exp = None
elif exp < 0:
exp = None
else:
exp * 7 // 365
if exp > exp_max:
exp = exp_max
return exp
# Marital status
_MariStat_dict = {'Other':0, 'Alone':1}
def transform_MariStat(self):
self._df['MariStat'] = self._df['MariStat'].map(self._MariStat_dict)
# Social category
def transform_SocioCateg(self):
self._df['SocioCateg'] = self._df['SocioCateg'].str.slice(0,4)
Подгружаем данные
data = InsDataFrame_Fr()
data.load_pd(df)
data.info()
Преобразовываем параметры
# Переименовываем
data.columns_match({'DrivAge':'driver_minage','LicAge':'driver_minexp'})
# Преобразовываем
data.transform_age()
data.transform_exp()
data.transform_gender()
data.transform_MariStat()
data.transform_SocioCateg()
# Пересечение пола и возраста, их квадраты
data.transform_age_gender()
data.polynomizer('driver_minage_m')
data.polynomizer('driver_minage_f')
Для переменных, содержащих более 2 значений:
# Onehot encoding
data.get_dummies(['VehUsage','SocioCateg'])
data.info()
col_features = [
'driver_minexp',
'Gender',
'MariStat',
'HasKmLimit',
'BonusMalus',
'OutUseNb',
'RiskArea',
'driver_minage_m',
'driver_minage_f',
'driver_minage_m_2',
'driver_minage_f_2',
'VehUsage_Private',
'VehUsage_Private+trip to office',
'VehUsage_Professional',
'VehUsage_Professional run',
'SocioCateg_CSP1',
'SocioCateg_CSP2',
'SocioCateg_CSP3',
'SocioCateg_CSP4',
'SocioCateg_CSP5',
'SocioCateg_CSP6',
'SocioCateg_CSP7'
]
col_target = ['ClaimAmount', 'ClaimsCount', 'AvgClaim']
df_freq = data.get_pd(col_features+col_target)
df_ac = df_freq[df_freq['ClaimsCount'] > 0].reset_index().copy()
from sklearn.model_selection import train_test_split
# Разбиение датасета для частоты на train/val/test
x_train_c, x_test_c, y_train_c, y_test_c = train_test_split(df_freq[col_features], df_freq.ClaimsCount, test_size=0.3, random_state=1)
x_valid_c, x_test_c, y_valid_c, y_test_c = train_test_split(x_test_c, y_test_c, test_size=0.5, random_state=1)
# Разбиение датасета для среднего убытка на train/val/test
x_train_ac, x_test_ac, y_train_ac, y_test_ac = train_test_split(df_ac[col_features], df_ac.AvgClaim, test_size=0.3, random_state=1)
x_valid_ac, x_test_ac, y_valid_ac, y_test_ac = train_test_split(x_test_ac, y_test_ac, test_size=0.5, random_state=1)
Градиентный бустинг - ансамблевый метод машинного обучения, использующийся для задач классификации, регрессии и ранжирования. Ансамбль представляет собой композицию простых базовых алгоритмов, в качестве которых обычно используются деревья решений.
Классический алгоритм GBM был предложен Джеромом Фридманом в 1999 году. Популярность методов GBM пришла в 2015-2016 гг. благодаря большому успеху библиотеки XGBoost в соревнованиях Kaggle.
Популярные библиотеки для GBM:
Все перечисленные библиотеки зачастую имеют сравнимый результат.
Модель
Пусть $y_i$ – значение переменной, которое необходимо предсказать, $x_i$ – входные данные.
Модель имеет вид $$\hat{y}_i = \sum_{k=1}^K f_k(x_i),\hspace{10pt} f_k \in \mathcal{F},$$ где $K$ – количество деревьев, $f$ – функция на пространстве $\mathcal{F}$, которое содержит все возможные деревья решений.
Целевая функция $$\text{Obj}(\theta) = L(\theta) + \Omega(\theta),$$ где
Тогда для представленной модели $\theta = \{f_1,f_2,\cdots,f_K\}$,
$$\text{Obj} = \sum_{i=1}^n l(y_i, \hat{y}_i) + \sum_{k=1}^K\Omega(f_k),$$где $l(y_i, \hat{y}_i)$ – функция потерь.
Обучение
Необходимо обучить функции $f_i$, каждая из которых включает структуру дерева и значения листьев.
Обозначим $\hat{y}_i^{(t)}$ предсказанное значение на шаге $t$. Тогда целевая функция имеет вид: $$\text{Obj} = \sum_{i=1}^n l\left(y_i, \hat{y}_i^{(t)}\right) + \sum_{i=1}^t\Omega(f_i).$$
Обучение деревьев происходит поочередно, начиная с постоянного предсказания: $$\begin{split}\hat{y}_i^{(0)} &= 0\\ \hat{y}_i^{(1)} &= f_1(x_i) = \hat{y}_i^{(0)} + f_1(x_i)\\ \hat{y}_i^{(2)} &= f_1(x_i) + f_2(x_i)= \hat{y}_i^{(1)} + f_2(x_i)\\ &\dots\\ \hat{y}_i^{(t)} &= \sum_{k=1}^t f_k(x_i)= \hat{y}_i^{(t-1)} + f_t(x_i)\end{split}$$
На каждом шаге выбирается дерево, которое оптимизирует целевую функцию. $$\begin{split}\text{Obj}^{(t)} & = \sum_{i=1}^n l\left(y_i, \hat{y}_i^{(t)}\right) + \sum_{i=1}^t\Omega(f_i) \\ & = \sum_{i=1}^n l\left(y_i, \hat{y}_i^{(t-1)} + f_t(x_i)\right) + \Omega(f_t) + \mathrm{constant}\end{split}$$
Для упрощения задачи оптимизации для заданной функции потерь используется разложение Тейлора: $$F(x+\Delta x) \simeq F(x) + F'(x)\Delta x + \frac{1}{2} F''(x)\Delta x^2 + \dots$$
Тогда обозначив градиент и гессиан функции потерь соответственно $$g_i = \partial_{\hat{y}_i^{(t-1)}} l(y_i, \hat{y}_i^{(t-1)}),\hspace{10pt}h_i = \partial_{\hat{y}_i^{(t-1)}}^2 l(y_i, \hat{y}_i^{(t-1)}),$$ целевая функция будет иметь вид $$\text{Obj}^{(t)} = \sum_{i=1}^n [l(y_i, \hat{y}_i^{(t-1)}) + g_i f_t(x_i) + \frac{1}{2} h_i f_t^2(x_i)] + \Omega(f_t) + \mathrm{constant}.$$
Убирая константы на шаге $t$, целевая функция упрощается в виде: $$\text{Obj}^{(t)} = \sum_{i=1}^n [g_i f_t(x_i) + \frac{1}{2} h_i f_t^2(x_i)] + \Omega(f_t).$$
Таким образом, величина потерь $L$ зависит только от $g_i$ и $h_i$.
Благодаря этому, XGBoost поддерживает пользовательские целевые функции, для которых достаточно задать градиент и гессиан функции потерь.
Регуляризация
Для начала определим дерево $f_t(x)$: $$f_t(x) = w_{q(x)}, w \in R^T, q:R^d\rightarrow \{1,2,\cdots,T\},$$ где $w$ – вектор значений на листьях дерева, $T$ – количество листьев, $q$ – функция, которая каждой точке набора данных ставит в соответствие лист дерева.
Тогда сложность модели имеет вид $$\Omega(f) = \gamma T + \frac{1}{2}\lambda \sum_{j=1}^T w_j^2+ \alpha \sum_{j=1}^T |w_j|,$$ где $\gamma$ – штраф на сложность деревьев, $\lambda$ – сила регуляризации $\ell_2$, $\alpha$ – сила регуляризации $\ell_1$.
Переобучение Для контроля переобучения помимо параметров $\gamma$, $\alpha$ и $\lambda$ используются также параметры:
max_depth)min_child_weight)subsample)colsample_bytree)Теория, стоящая за оценкой весов на листьях и нахождением разделений выходит за рамки нашего рассмотрения.
Основные гиперпараметры, используемые в библиотеке XGBoost
Стратегии подбора гиперпараметров:
Использование некоторых алгоритмов подбора гиперпараметров для XGBoost
Для подбора параметров воспользуемся реализацией алгоритма Tree-structured Parzen Estimator (TPE) в библиотеке Hyperopt. Алгоритм использует подход последовательной оптимизации, основанной на модели (sequential model-based optimization, SMBO). Метод основывается в байесовской оптимизации и гауссовских процессах.
#!pip install hyperopt --upgrade
from functools import partial
#Anaconda Prompt command:
#conda install xgboost -- does not found
#pip install xgboost -- success
import xgboost as xgb
from hyperopt import hp, fmin, tpe, space_eval, Trials, STATUS_OK
# Конвертация наборов данных в формат, поддерживающийся XGBoost
train_c = xgb.DMatrix(x_train_c, y_train_c)
valid_c = xgb.DMatrix(x_valid_c, y_valid_c)
test_c = xgb.DMatrix(x_test_c, y_test_c)
# Зададим функцию Deviance для распределения Пуассона
def xgb_eval_dev_poisson(yhat, y):
t_hat, t = yhat + 1, y.get_label() + 1
return 'dev_poisson', 2 * np.sum(t * np.log(t / t_hat) - (t - t_hat))
# Определим функцию для оптимизации гиперпараметров алгоритмом TPE
def objective(params, cv_params, data):
if 'max_depth' in params.keys():
params['max_depth'] = int(params['max_depth'])
cv_result = xgb.cv(params=params, dtrain=data, **cv_params)
name = [i for i in cv_result.columns if all([i.startswith('test-'), i.endswith('-mean')])][-1]
score = cv_result[name][-1:].values[0]
return {'loss': score, 'status': STATUS_OK}
# Определим параметры выполнения кроссвалидации
cv_params = {'num_boost_round': 300,
'nfold': 5,
'shuffle': True,
'stratified': False,
'feval': xgb_eval_dev_poisson,
'maximize': False,
'early_stopping_rounds': 20
}
# Определим границы, в которых будем искать гиперпараметры
space_freq = {'objective': 'count:poisson',
'max_depth': hp.choice('max_depth', [5, 8, 10, 12, 15]),
'min_child_weight': hp.uniform('min_child_weight', 0, 50),
'subsample': hp.uniform('subsample', 0.5, 1),
'colsample_bytree': hp.uniform('colsample_bytree', 0.5, 1),
'alpha': hp.uniform('alpha', 0, 1),
'lambda': hp.uniform('lambda', 0, 1),
'eta': hp.uniform('eta', 0.01, 1),
'tree_method': 'hist'
}
# Оптимизация (количество итераций снижено для ускорения работы)
trials = Trials()
best = fmin(fn=partial(objective, cv_params=cv_params, data=train_c),
space=space_freq, trials=trials, algo=tpe.suggest, max_evals=50, timeout=3600)
# Оптимальные гиперпараметры
best_params = space_eval(space_freq, best)
best_params
train_params = {'num_boost_round': 300,
'feval': xgb_eval_dev_poisson,
'maximize': False,
'verbose_eval': False}
# Построение модели с ранней остановкой (early stopping)
progress = dict()
xgb_freq = xgb.train(params=best_params, dtrain=train_c, early_stopping_rounds=10,
evals=[(train_c, "train"), (valid_c, "valid")],
evals_result=progress, **train_params)
# Построение модели без ранней остановки
progress_wo_es = dict()
xgb_freq_wo_es = xgb.train(params=best_params, dtrain=train_c, evals=[(train_c, "train"), (valid_c, "valid"), (test_c, "test")],
evals_result=progress_wo_es, **train_params)
import matplotlib.pyplot as plt
plt.subplots(2,2, figsize=(10,6))
plt.subplot(2,2,1)
plt.plot(progress_wo_es['valid']['poisson-nloglik'], label='valid')
plt.plot(progress['valid']['poisson-nloglik'], label='valid early stopping', linestyle='dashed', color='black')
plt.plot(progress_wo_es['test']['poisson-nloglik'], label='test')
plt.xlabel('Epochs'); plt.ylabel('Negative Log-Likelihood'); plt.legend(); plt.tight_layout()
plt.subplot(2,2,2)
plt.plot(progress_wo_es['valid']['dev_poisson'], label='valid')
plt.plot(progress['valid']['dev_poisson'], label='valid early stopping', linestyle='dashed', color='black')
plt.plot(progress_wo_es['test']['dev_poisson'], label='test')
plt.xlabel('Epochs'); plt.ylabel('Poisson Deviance'); plt.legend(); plt.tight_layout()
plt.subplot(2,2,3)
plt.plot(progress_wo_es['valid']['poisson-nloglik'], label='valid')
plt.plot(progress['valid']['poisson-nloglik'], label='valid early stopping', linestyle='dashed', color='black')
plt.plot(progress_wo_es['test']['poisson-nloglik'], label='test')
plt.plot(progress_wo_es['train']['poisson-nloglik'], label='train')
plt.xlabel('Epochs'); plt.ylabel('Negative Log-Likelihood'); plt.legend(); plt.tight_layout()
plt.subplot(2,2,4)
plt.plot(progress_wo_es['valid']['dev_poisson'], label='valid')
plt.plot(progress['valid']['dev_poisson'], label='valid early stopping', linestyle='dashed', color='black')
plt.plot(progress_wo_es['test']['dev_poisson'], label='test')
plt.plot(progress_wo_es['train']['dev_poisson'], label='train')
plt.xlabel('Epochs'); plt.ylabel('Poisson Deviance'); plt.legend(); plt.tight_layout(); plt.show()
# Отбор признаков (Feature Importance)
importance_type = ['total_gain', 'gain', 'weight', 'total_cover', 'cover']
xgb.plot_importance(xgb_freq, importance_type=importance_type[1]); plt.show()
# Конвертация наборов данных в формат, поддерживающийся XGBoost
train_ac = xgb.DMatrix(x_train_ac, y_train_ac)
valid_ac = xgb.DMatrix(x_valid_ac, y_valid_ac)
test_ac = xgb.DMatrix(x_test_ac, y_test_ac)
# Зададим функцию Deviance для гамма-распределения
def xgb_eval_dev_gamma(yhat, dtrain):
y = dtrain.get_label()
return 'dev_gamma', 2 * np.sum(-np.log(y/yhat) + (y-yhat)/yhat)
# Определим границы, в которых будем искать гиперпараметры
space_avgclm = {'objective': 'reg:gamma',
'max_depth': hp.choice('max_depth', [5, 8, 10, 12, 15]),
'min_child_weight': hp.uniform('min_child_weight', 0, 50),
'subsample': hp.uniform('subsample', 0.5, 1),
'colsample_bytree': hp.uniform('colsample_bytree', 0.5, 1),
'alpha': hp.uniform('alpha', 0, 1),
'lambda': hp.uniform('lambda', 0, 1),
'eta': hp.uniform('eta', 0.01, 1),
'tree_method': 'hist'
}
# Определим параметры выполнения кроссвалидации
cv_params = {'num_boost_round': 300,
'nfold': 5,
'shuffle': True,
'stratified': False,
'feval': xgb_eval_dev_gamma,
'maximize': False,
'early_stopping_rounds': 20
}
# Оптимизация (количество итераций снижено для ускорения работы)
trials = Trials()
best = fmin(fn=partial(objective, cv_params=cv_params, data=train_ac),
space=space_avgclm, trials=trials, algo=tpe.suggest, max_evals=50, timeout=3600)
# Оптимальные гиперпараметры
best_params = space_eval(space_avgclm, best)
best_params
train_params = {'num_boost_round': 300,
'feval': xgb_eval_dev_gamma,
'maximize': False,
'verbose_eval': False}
# Построение модели с ранней остановкой (early stopping)
progress = dict()
xgb_avgclaim = xgb.train(params=best_params, dtrain=train_ac, early_stopping_rounds=10, evals=[(train_ac, "train"), (valid_ac, "valid")],
evals_result=progress, **train_params)
# Построение модели без ранней остановки
progress_wo_es = dict()
xgb_avgclaim_wo_es = xgb.train(params=best_params, dtrain=train_ac, evals=[(train_ac, "train"), (valid_ac, "valid"), (test_ac, "test")],
evals_result=progress_wo_es, **train_params)
plt.subplots(2,2, figsize=(10,6))
plt.subplot(2,2,1)
plt.plot(progress_wo_es['valid']['gamma-nloglik'], label='valid')
plt.plot(progress['valid']['gamma-nloglik'], label='valid early stopping', linestyle='dashed', color='black')
plt.plot(progress_wo_es['test']['gamma-nloglik'], label='test')
plt.xlabel('Epochs'); plt.ylabel('Negative Log-Likelihood'); plt.legend(); plt.tight_layout()
plt.subplot(2,2,2)
plt.plot(progress_wo_es['valid']['dev_gamma'], label='valid')
plt.plot(progress['valid']['dev_gamma'], label='valid early stopping', linestyle='dashed', color='black')
plt.plot(progress_wo_es['test']['dev_gamma'], label='test')
plt.xlabel('Epochs'); plt.ylabel('Gamma Deviance'); plt.legend(); plt.tight_layout()
plt.subplot(2,2,3)
plt.plot(progress_wo_es['valid']['gamma-nloglik'], label='valid')
plt.plot(progress['valid']['gamma-nloglik'], label='valid early stopping', linestyle='dashed', color='black')
plt.plot(progress_wo_es['test']['gamma-nloglik'], label='test')
plt.plot(progress_wo_es['train']['gamma-nloglik'], label='train')
plt.xlabel('Epochs'); plt.ylabel('Negative Log-Likelihood'); plt.legend(); plt.tight_layout()
plt.subplot(2,2,4)
plt.plot(progress_wo_es['valid']['dev_gamma'], label='valid')
plt.plot(progress['valid']['dev_gamma'], label='valid early stopping', linestyle='dashed', color='black')
plt.plot(progress_wo_es['test']['dev_gamma'], label='test')
plt.plot(progress_wo_es['train']['dev_gamma'], label='train')
plt.xlabel('Epochs'); plt.ylabel('Gamma Deviance'); plt.legend(); plt.tight_layout(); plt.show()
# Отбор признаков (Feature Importance)
importance_type = ['total_gain', 'gain', 'weight', 'total_cover', 'cover']
xgb.plot_importance(xgb_avgclaim, importance_type=importance_type[1]); plt.show()
#!pip install eli5
import eli5
eli5.show_weights(xgb_avgclaim)
eli5.explain_prediction(xgb_avgclaim, x_test_ac.iloc[0,:])
np.log(xgb_avgclaim.predict(xgb.DMatrix(x_test_ac.iloc[[0],:])))
#!pip install shap
#Anaconda Prompt:
#conda install shap -- does not install
#pip install shap -- does not install
#conda install -c conda-forge shap --success!
#conda-forge is a community effort that provides conda packages for a wide range of software.
#Missing a package that you would love to install with conda? - Chances are we have already packaged it for you!
#pip install --upgrade setuptools
import shap
explainer = shap.TreeExplainer(xgb_avgclaim)
shap_values = explainer.shap_values(x_test_ac)
shap.initjs()
shap.force_plot(explainer.expected_value, shap_values[0,:], x_test_ac.iloc[0,:])
shap.waterfall_plot(explainer.expected_value, shap_values[0,:], x_test_ac.iloc[0,:])
shap.initjs()
shap.force_plot(explainer.expected_value, shap_values, x_test_ac)
shap.dependence_plot("BonusMalus", shap_values, x_test_ac)
shap.summary_plot(shap_values, x_test_ac)
shap.summary_plot(shap_values, x_test_ac, plot_type="bar")
predictions = pd.DataFrame()
predictions['CountPredicted'] = xgb_freq.predict(xgb.DMatrix(df_freq[col_features]))
predictions['AvgClaimPredicted'] = xgb_avgclaim.predict(xgb.DMatrix(df_freq[col_features]))
predictions['CountPredicted'].min()
predictions['BurningCost'] = predictions.CountPredicted * predictions.AvgClaimPredicted
predictions.head()
Об особенностях сохранения моделей:
xgb_avgclaim.save_model('D:\Cloud\Git\geekbrains-ml-business\\avg_claim.model')
xgb_avgclaim.feature_names
xgb_avgclaim.feature_types
model =xgb.Booster()
model.load_model('D:\Cloud\Git\geekbrains-ml-business\\avg_claim.model')
print(model.feature_names)
type(model.feature_names)
model.feature_types
В текущем домашнем задание предлагается взглянуть на задачу моделирования количества страховых случаев как на задачу многоклассовой классификации.
#pip install "git+https://github.com/MindSetLib/Insuranse_Lesson.git#egg=insolver"
import pandas as pd
import numpy as np
import xgboost as xgb
import matplotlib.pyplot as plt
from sklearn.metrics import f1_score, accuracy_score, confusion_matrix
#from insolver import InsDataframe, InsolverGradientBoostingWrapper, train_val_test_split
import warnings
from functools import partial
from xgboost import DMatrix, cv as xcv, train as xtrain, XGBClassifier, XGBRegressor
#from lightgbm import Dataset, cv as lcv, train as ltrain, LGBMClassifier, LGBMRegressor
#from catboost import Pool, cv as ccv, train as ctrain, CatBoostClassifier, CatBoostRegressor
from hyperopt import hp, tpe, space_eval, Trials, STATUS_OK, STATUS_FAIL, fmin
from sklearn.model_selection import train_test_split
#взято отсюда https://github.com/MindSetLib/Insuranse_Lesson/blob/master/insolver/core.py
def train_val_test_split(x, y, val_size, test_size, random_state=0, shuffle=True, stratify=None):
x_train, x_test, y_train, y_test = train_test_split(x, y, random_state=random_state, shuffle=shuffle,
test_size=test_size, stratify=stratify)
x_train, x_valid, y_train, y_valid = train_test_split(x_train, y_train, random_state=random_state, shuffle=shuffle,
test_size=val_size/(1-test_size), stratify=stratify)
return x_train, x_valid, x_test, y_train, y_valid, y_test
def objective_gb(params, algorithm, cv_params, data_params, X, y):
if data_params is None:
data_params = dict()
for x in ['max_depth', 'num_boost_round', 'max_leaves', 'max_bin', 'num_leaves',
'min_child_samples', 'min_data_in_leaf']:
if x in params.keys():
params[x] = int(params[x])
if algorithm == 'xgboost':
dtrain = DMatrix(X, y, **data_params)
cv_result = xcv(params=params, dtrain=dtrain, **cv_params)
name = [i for i in cv_result.columns if all([i.startswith('test-'), i.endswith('-mean')])][-1]
score = cv_result[name][-1:].values[0]
elif algorithm == 'lightgbm':
dtrain = Dataset(X, y, **data_params)
cv_result = lcv(params=params, train_set=dtrain, **cv_params)
name = [i for i in cv_result.keys() if i.endswith('-mean')][-1]
score = cv_result[name][-1]
elif algorithm == 'catboost':
dtrain = Pool(X, y, **data_params)
cv_result = ccv(params=params, dtrain=dtrain, **cv_params)
name = [i for i in cv_result.columns if all([i.startswith('test-'), i.endswith('-mean')])][-1]
score = cv_result[name][-1:].values[0]
else:
warnings.warn('Error occurred in "algorithm" attribute')
score = 0
return {'loss': score, 'status': STATUS_OK}
#взято отсюда: https://github.com/MindSetLib/Insuranse_Lesson/blob/master/insolver/core.py
class InsolverGradientBoostingWrapper(object):
def __init__(self, algorithm, task='classification'):
if algorithm in ['xgboost', 'lightgbm', 'catboost']:
self.algorithm = algorithm
else:
warnings.warn('Specified algorithm parameter is not supported. '
'Try to enter one of the following options: ["xgboost", "lightgbm", "catboost"].')
if task in ['regression', 'classification']:
self.task = task
else:
warnings.warn('Specified task parameter is not supported. '
'Try to enter one of the following options: ["regression", "classification"].')
self.trials, self.best_params, self.core_params, self.data_params = None, None, None, None
self.model, self.booster = None, None
@staticmethod
def cv_parameters_default_xgboost():
return {'num_boost_round': 10,
'nfold': 3,
'stratified': False,
'folds': None,
'metrics': (),
'obj': None,
'feval': None,
'maximize': False,
'early_stopping_rounds': None,
'fpreproc': None,
'as_pandas': True,
'verbose_eval': None,
'show_stdv': True,
'seed': 0,
'callbacks': None,
'shuffle': True}
@staticmethod
def cv_parameters_default_lightgbm():
return {'num_boost_round': 100,
'folds': None,
'nfold': 5,
'stratified': True,
'shuffle': True,
'metrics': None,
'fobj': None,
'feval': None,
'init_model': None,
'feature_name': 'auto',
'categorical_feature': 'auto',
'early_stopping_rounds': None,
'fpreproc': None,
'verbose_eval': None,
'show_stdv': True,
'seed': 0,
'callbacks': None,
'eval_train_metric': False,
'return_cvbooster': False}
@staticmethod
def cv_parameters_default_catboost():
return {'iterations': None,
'num_boost_round': None,
'fold_count': 3,
'nfold': None,
'inverted': False,
'partition_random_seed': 0,
'seed': None,
'shuffle': True,
'logging_level': None,
'stratified': None,
'as_pandas': True,
'metric_period': None,
'verbose': None,
'verbose_eval': None,
'plot': False,
'early_stopping_rounds': None,
'folds': None,
'type': 'Classical'}
def hyperopt_cv(self, X, y, params, cv_params, data_params=None, max_evals=10, fn=None, algo=None, timeout=None):
trials = Trials()
if data_params is not None:
self.data_params = data_params
if algo is None:
algo = tpe.suggest
if fn is None:
fn = partial(objective_gb, algorithm=self.algorithm, X=X, y=y, cv_params=cv_params, data_params=data_params)
try:
best = fmin(fn=fn, space=params, trials=trials, algo=algo, max_evals=max_evals, timeout=timeout)
best_params = space_eval(params, best)
for x in ['max_depth', 'num_boost_round', 'max_leaves', 'max_bin', 'num_leaves',
'min_child_samples', 'min_data_in_leaf']:
if x in best_params.keys():
best_params[x] = int(best_params[x])
self.best_params, self.trials = best_params, trials
core_params = ['num_boost_round', 'obj', 'feval', 'maximize', 'early_stopping_rounds',
'verbose_eval', 'callbacks', 'fobj', 'init_model', 'feature_name', 'categorical_feature',
'iterations', 'logging_level', 'metric_period', 'verbose', 'plot']
self.core_params = {i: cv_params[i] for i in cv_params if i in core_params}
except Exception as e:
return {'status': STATUS_FAIL, 'exception': str(e)}
def fit_booster(self, X, y, data_params=None, core_params=None):
dtrain_params, train_params = dict(), dict()
if self.data_params is not None:
dtrain_params.update(self.data_params)
if data_params is not None:
dtrain_params.update(data_params)
if self.core_params is not None:
train_params.update(self.core_params)
if core_params is not None:
train_params.update(core_params)
if self.best_params is not None:
params = self.best_params
else:
params = {}
if self.algorithm == 'xgboost':
dtrain = DMatrix(X, y, **dtrain_params)
if 'evals' in train_params.keys():
train_params['evals'] = [(DMatrix(i[0][0], i[0][1]), i[1]) for i in train_params['evals']]
self.booster = xtrain(params=params, dtrain=dtrain, **train_params)
elif self.algorithm == 'lightgbm':
dtrain = Dataset(X, y, **dtrain_params)
if 'evals' in train_params.keys():
evals = train_params.pop('evals')
train_params['valid_names'] = [i[1] for i in evals]
train_params['valid_sets'] = [(Dataset(i[0][0], i[0][1])) for i in evals]
self.booster = ltrain(params=params, train_set=dtrain, **train_params)
elif self.algorithm == 'catboost':
dtrain = Pool(X, y, **dtrain_params)
if 'evals' in train_params.keys():
train_params['evals'] = [Pool(i[0][0], i[0][1]) for i in train_params['evals']]
self.booster = ctrain(params=params, dtrain=dtrain, **train_params)
else:
warnings.warn('Specified algorithm parameter is not supported.')
def model_init(self, params=None):
hparams = dict()
if self.core_params is not None:
hparams.update(self.core_params)
if self.best_params is not None:
hparams.update(self.best_params)
if params is not None:
hparams.update(params)
aliases = {'eta': 'learning_rate', 'boosting': 'boosting_type', 'max_leaves': 'num_leaves',
'num_iterations': 'n_estimators', 'num_iteration': 'n_estimators',
'n_iter': 'n_estimators', 'num_tree': 'n_estimators', 'num_trees': 'n_estimators',
'num_round': 'n_estimators', 'num_rounds': 'n_estimators', 'num_boost_round': 'n_estimators'}
for x in hparams.keys():
if x in aliases.keys():
hparams[aliases[x]] = hparams.pop(x)
if self.task == 'classification':
if self.algorithm == 'xgboost':
self.model = XGBClassifier(**hparams)
elif self.algorithm == 'lightgbm':
self.model = LGBMClassifier(**hparams)
elif self.algorithm == 'catboost':
self.model = CatBoostClassifier(**hparams)
else:
warnings.warn('Specified algorithm parameter is not supported.')
elif self.task == 'regression':
if self.algorithm == 'xgboost':
self.model = XGBRegressor(**hparams)
elif self.algorithm == 'lightgbm':
self.model = LGBMRegressor(**hparams)
elif self.algorithm == 'catboost':
self.model = CatBoostRegressor(**hparams)
else:
warnings.warn('Specified algorithm parameter is not supported.')
else:
warnings.warn('Tasks other than "classification" and "regression" are not supported.')
def fit(self, X, y, **kwargs):
if self.model is None:
warnings.warn('Model is not initiated, please use .model_init() method.')
else:
self.model.fit(X, y, **kwargs)
#df = pd.read_csv('/content/drive/My Drive/Colab Notebooks/freMPL-R.csv', low_memory=False)
df = pd.read_csv('D:\Cloud\Git\geekbrains-ml-business\\freMPL-R.csv', low_memory=False)
df = df.loc[df.Dataset.isin([5, 6, 7, 8, 9])]
df.drop('Dataset', axis=1, inplace=True)
df.dropna(axis=1, how='all', inplace=True)
df.drop_duplicates(inplace=True)
df.reset_index(drop=True, inplace=True)
df.info()
Предобработайте данные
В предыдущем уроке мы заметили отрицательную величину убытка для некоторых наблюдений. Заметим, что для всех таких полисов переменная "ClaimInd" принимает только значение 0. Поэтому заменим все соответствующие значения "ClaimAmount" нулями.
NegClaimAmount = df.loc[df.ClaimAmount < 0, ['ClaimAmount','ClaimInd']]
print('Unique values of ClaimInd:', NegClaimAmount.ClaimInd.unique())
df.loc[df.ClaimAmount < 0, 'ClaimAmount'] = 0
df.info()
df.head()
Для моделирования частоты убытков сгенерируем показатель как сумму индикатора того, что убыток произошел ("ClaimInd") и количества заявленных убытков по различным видам ущерба за 4 предшествующих года ("ClaimNbResp", "ClaimNbNonResp", "ClaimNbParking", "ClaimNbFireTheft", "ClaimNbWindscreen").
В случаях, если соответствующая величина убытка равняется нулю, сгенерированную частоту также обнулим.
df['ClaimsCount'] = df.ClaimInd + df.ClaimNbResp + df.ClaimNbNonResp + df.ClaimNbParking + df.ClaimNbFireTheft + df.ClaimNbWindscreen
df.loc[df.ClaimAmount == 0, 'ClaimsCount'] = 0
df.drop(["ClaimNbResp", "ClaimNbNonResp", "ClaimNbParking", "ClaimNbFireTheft", "ClaimNbWindscreen"], axis=1, inplace=True)
pd.DataFrame(df.ClaimsCount.value_counts()).rename({'ClaimsCount': 'Policies'}, axis=1)
data = InsDataFrame_Fr()
data.load_pd(df)
data.info()
# Переименовываем
data.columns_match({'DrivAge':'driver_minage','LicAge':'driver_minexp'})
# Преобразовываем
data.transform_age()
data.transform_exp()
data.transform_gender()
data.transform_MariStat()
data.transform_SocioCateg()
# Пересечение пола и возраста, их квадраты
data.transform_age_gender()
data.polynomizer('driver_minage_m')
data.polynomizer('driver_minage_f')
# Onehot encoding Для переменных, содержащих более 2 значений:
data.get_dummies(['VehUsage','SocioCateg'])
data.info()
col_features = [
'driver_minexp',
'Gender',
'MariStat',
'HasKmLimit',
'BonusMalus',
'OutUseNb',
'RiskArea',
'driver_minage_m',
'driver_minage_f',
'driver_minage_m_2',
'driver_minage_f_2',
'VehUsage_Private',
'VehUsage_Private+trip to office',
'VehUsage_Professional',
'VehUsage_Professional run',
'SocioCateg_CSP1',
'SocioCateg_CSP2',
'SocioCateg_CSP3',
'SocioCateg_CSP4',
'SocioCateg_CSP5',
'SocioCateg_CSP6',
'SocioCateg_CSP7'
]
#col_target = ['ClaimAmount', 'ClaimsCount', 'AvgClaim']
col_target = ['ClaimAmount','ClaimsCount']
df_freq = data.get_pd(col_features+col_target)
df_freq.info()
XGBoost для многоклассовой классификации принимает на вход значения меток классов в виде [0, num_classes]. Заменим значение 11 на 10.
df_freq['ClaimsCount'] = df_freq['ClaimsCount'].replace({11: 10}).fillna(0)
Посмотрим, сколько полисов соответствуют каждому из значений ClaimsCount, используя метод groupby. Для полученных значений также посчитаем нормированную частоту.
FreqCount = pd.DataFrame(df_freq.groupby('ClaimsCount').size(), columns=['Count'])
#FreqCount['Freq'] = '<Ваш код здесь>'
FreqCount.Count.plot(kind='bar')
plt.ylabel('Frequency')
plt.show()
FreqCount
Заметим, что в данном случае присутствует проблема несбалансированности классов. Поэтому, для того, чтобы по возможности избежать ее, воспользуемся взвешиванием наблюдений для обучения модели. Для этого в исходном наборе данных создадим столбец weight. Присвоим ему некоторые значения, например, можно задать 0.05 для значений ClaimsCount 0, а для остальных - 1 (Для этого можем использовать функцию np.where). Также можно попробовать какой-либо другой способ задания весов, приведенный пример не гарантирует хороших результатов.
#поставим свое значение
df_freq['weight'] = np.where(df_freq['ClaimsCount'] == 0, 0.15, 1)
df_freq
Разобьем имеющийся набор данных на обучающую, валидационную и тестовую выборки в отношениях 70%/15%/15% соответственно. Зададим зерно для случайного разбиения равным 10.
x_train, x_valid, x_test, \
y_train, y_valid, y_test = train_val_test_split(df_freq.drop(['ClaimAmount','ClaimsCount'], axis = 1)
,df_freq.ClaimsCount,val_size=0.15, test_size=0.15, random_state=10)
Далее, создадим объекты DMatrix для обучающей, валидационной и тестовой выборок. Для обучающей выборки также укажем параметр weight равным полученному ранее столбцу весов. Данный столбец также нужно исключить из объекта передаваемого в параметр data.
igb = InsolverGradientBoostingWrapper(algorithm='xgboost')
space_xgboost = {'objective': 'multi:softmax',
'max_depth': hp.choice('max_depth', [5, 8, 10, 12, 15]),
'min_child_weight': hp.uniform('min_child_weight', 0, 50),
'subsample': hp.uniform('subsample', 0.5, 1),
'colsample_bytree': hp.uniform('colsample_bytree', 0.5, 1),
'alpha': hp.uniform('alpha', 0, 1),
'lambda': hp.uniform('lambda', 0, 1),
'eta': hp.uniform('eta', 0.01, 1),
'gamma': hp.uniform('gamma', 0.01, 1000),
'num_class': len(df_freq.ClaimsCount.unique()),#тут число классов
'tree_method': 'hist'
}
cv_params = {'num_boost_round': 1000,
'nfold': 3,
'early_stopping_rounds': 20,
'seed': 0,
'shuffle': True,
'stratified': False }
igb.hyperopt_cv(x_train.drop('weight', axis=1), y_train, space_xgboost, cv_params,
data_params={'weight': x_train['weight']}, max_evals=50)
Для оптимизации гиперпараметров можно воспользоваться различными методами.
Далее обучим нашу модель с оптимальными параметрами
igb.fit_booster(x_train.drop('weight', axis=1), y_train, core_params={'evals': [((x_train.drop('weight', axis=1),y_train),'train'),
((x_valid.drop('weight', axis=1),y_valid),'valid')]})
Посчитаем метрики accuracy и f1 на наших наборах данных, также можем визуализировать confusion matrix, например, с помощью plt.imshow(). Можно использовать предложенный ниже код.
dfsets = [{'set': 'train', 'x': x_train.drop('weight', axis=1), 'target': y_train},
{'set': 'valid', 'x': x_valid.drop('weight', axis=1), 'target': y_valid},
{'set': 'test', 'x': x_test.drop('weight', axis=1), 'target': y_test}]
for dfset in dfsets:
class_preds = igb.booster.predict(xgb.DMatrix(dfset['x'])) # Посчитаем предсказанные значения
print('F1 Score on ' + str(dfset['set'])+':', f1_score(dfset['target'],class_preds, average='micro')) # Посчитаем F1 Score
plt.subplots(1,3, figsize=(15,3))
for i in range(len(dfsets)):
confmatrix = confusion_matrix(dfsets[i]['target'], igb.booster.predict(xgb.DMatrix(dfsets[i]['x'])))
plt.subplot(1,3,i+1)
plt.imshow(confmatrix, cmap='Greys')
plt.colorbar()
plt.ylabel('True')
plt.xlabel('Predicted')
plt.show()
Как вы оцениваете качество построенной модели? Какие проблемы могут здесь присутствовать? Как можно улучшить результат?
Ответ: изначально наблюдался дисбаланс классов (возобладал класс 0), отсюда такой результат. Подобную задачу следует решать с помощью регрессии, а не классификации.